# 加载相关包
library(openxlsx)
library(geepack)
library(gtsummary)
library(mice)
# 加载自定义函数
source("./func/gee_run.r")

############# 需求 #############
# 和 gee_918.R 基本一致，但是修改了部分操作
# 1. 多因素删除缺失值
#    - 分别剔除因变量和自变量中有缺失值的行 => 一次性剔除有缺失值的行（让结果更可靠）
# 2. 单因素虚拟填补缺失值
#    - 
# 3. 多因素 mice 填补缺失值
#    - 使用 mice 填补因变量和协变量的缺失值 => 删除因变量缺失的行，填补协变量的缺失值
#    - 输出信息：每个协变量填补前后各有多少行（观察填补行数是否过多导致结果不可靠）
############# 需求 #############

###### 读取数据，数据预处理
data_test = read.xlsx("./input/data918.xlsx", 3)
data_anal = data_test[, c(1:2, 9:13, 14:25, 27, 29:39)]
varsx = colnames(data_test)[c(14:25, 27)]  # 自变量
varsy = colnames(data_test)[c(9:13)]       # 因变量
varsc = colnames(data_test)[c(29:39)]      # 协变量

# 确保变量类型一致：因变量为数值，自变量、协变量为因子
data_anal$id = as.character(data_anal$id)
data_anal[varsy] = lapply(data_anal[varsy], as.numeric)
data_anal[varsx] = lapply(data_anal[varsx], as.factor)
data_anal[varsc] = lapply(data_anal[varsc], as.factor)
# str(data_anal)
######


# 1. 多因素删除缺失值
## 1.1 先删除包含缺失值的行 
tmp_data = data_anal[complete.cases(data_anal), ]
tmp_data[varsx] = lapply(tmp_data[varsx], as.numeric)  # 删除缺失值后，某些因子项的 level 计数可能变为 0，因此需要剔除该 level
tmp_data[varsc] = lapply(tmp_data[varsc], as.numeric)  # 同上
tmp_data[varsx] = lapply(tmp_data[varsx], as.factor)
tmp_data[varsc] = lapply(tmp_data[varsc], as.factor)
## 1.2
gee_run(data = tmp_data,
        varsx = varsx,
        varsy = varsy,
        varsc = varsc[-5],
        file = "./output/0919_delete_10cv_del_birthweight_c.xlsx")
gee_run(data = tmp_data,
        varsx = varsx,
        varsy = varsy,
        varsc = varsc[-2],
        file = "./output/0919_delete_10cv_del_week_c.xlsx")



# 2. 单因素虚拟填补缺失值
## 2.1 预处理，将缺失值定义为该列的最大值 + 1
tmp_data = data_anal
tmp_data[varsx] = lapply(tmp_data[varsx], as.numeric)  # 先转换为数值，才能取得到最大值
for (k in varsx) {
  tmp_data[, k] = ifelse(is.na(tmp_data[, k]), max(tmp_data[, k], na.rm = T), tmp_data[, k])
}
# tmp_data[is.na(tmp_data)] = 99                                    # 或者快捷设置为一个特定值
# tmp_data[is.na(tmp_data)] = max(temp_data[varsx], na.rm = T) + 1  # 或者快捷设置为数据中没有的值
tmp_data[varsx] = lapply(tmp_data[varsx], as.factor)   # 再转换回因子
## 2.2 gee（单因素不要传入 varsc 参数）
gee_run(tmp_data,
        varsx = varsx,
        varsy = varsy,
        file = "./output/0919_virtual.xlsx")



# 3. 多因素 mice 填补缺失值
## 3.1 先跑一次 mice
tmp_data = unique(data_anal[, c("id", varsc)])
get_na_unit = function (x) {
  t1 = data.frame(id = tmp_data$id, var = x)
  t2 = t1[complete.cases(t1), ]
  return (length(unique(t2$id)))
}
total = length(unique(tmp_data$id))
##### mice_before
no_na_unit = as.vector(unlist(lapply(tmp_data, get_na_unit)))
has_na_unit = as.vector(unlist(lapply(no_na_unit, function (x) total - x )))
mice_before = data.frame(variable = colnames(tmp_data), no_na_unit, has_na_unit)
mice_before$total = total
#####

miss = md.pattern(tmp_data)                # 查看缺失数据的模式
imp = mice(tmp_data, m = 5, seed = 1)      # 填补
mice_data = complete(imp, action = 1)      # 懒惰的人选择第一个
mice_data = merge(data_anal[, c("id", "month", varsy, varsx)], mice_data)   # 合并


##### 每个自变量的缺失个体
tmp_data = data_anal[, c("id", varsx)]
total = length(unique(tmp_data$id))
no_na_x_unit = as.vector(unlist(lapply(tmp_data, get_na_unit)))
has_na_x_unit = as.vector(unlist(lapply(no_na_x_unit, function (x) total - x )))
x_missing = data.frame(variable = colnames(tmp_data), no_na_x_unit, has_na_x_unit)
x_missing$total = total


##### 自变量不缺失的个体中，协变量缺失情况
tmp_data = data_anal[, c("id", varsx, varsc)]
tmp_data = tmp_data[complete.cases(tmp_data[, varsx]), ]
total = length(unique(tmp_data$id))
no_na_c_unit = as.vector(unlist(lapply(tmp_data[, c("id", varsc)], get_na_unit)))
has_na_c_unit = as.vector(unlist(lapply(no_na_c_unit, function (x) total - x )))
c_missing = data.frame(variable = colnames(tmp_data[, c("id", varsc)]), no_na_c_unit, has_na_c_unit)
c_missing$total = total

wb = createWorkbook()
addWorksheet(wb, sheet = "1372_x_missing")
writeDataTable(wb, sheet = "1372_x_missing", x = x_missing)       # 写入自变量的缺失情况
addWorksheet(wb, sheet = "767_c_missing")
writeDataTable(wb, sheet = "767_c_missing", x = c_missing)       # 写入自变量的缺失情况
addWorksheet(wb, sheet = "1372_c_mice_before")
writeDataTable(wb, sheet = "1372_c_mice_before", x = mice_before)       # 写入 mice 前后的 varsc 数据变化
addWorksheet(wb, sheet = "mice_data")
writeDataTable(wb, sheet = "mice_data", x = mice_data)           # 写入 mice 后的数据
saveWorkbook(wb, "./output/0919_mice_data.xlsx", overwrite = T)  # 输出到 xlsx

mice_data = read.xlsx("./output/0919_mice_data.xlsx", 4)
mice_data = mice_data[complete.cases(mice_data[, varsx]), ]
## 3.2 gee
gee_run(mice_data,
        varsx = varsx,
        varsy = varsy,
        varsc = varsc[-5],
        file = "./output/0919_mice_10cv_del_birthweight_c.xlsx")
gee_run(mice_data,
        varsx = varsx,
        varsy = varsy,
        varsc = varsc[-2],
        file = "./output/0919_mice_10cv_del_week_c.xlsx")
